Single-nucleus multiomics reveals the disrupted regulatory programs in three brain regions of sporadic early-onset Alzheimer’s disease

Sporadic early-onset Alzheimer’s disease (sEOAD) represents a significant but less-studied subtype of Alzheimer’s disease (AD). Here, we generated a single-nucleus multiome atlas derived from the postmortem prefrontal cortex, entorhinal cortex, and hippocampus of nine individuals with or without sEOAD. Comprehensive analyses were conducted to delineate cell type-specific transcriptomic changes and linked candidate cis-regulatory elements (cCREs) across brain regions. We prioritized seven conservative transcription factors in glial cells in multiple brain regions, including RFX4 in astrocytes and IKZF1 in microglia, which are implicated in regulating sEOAD-associated genes. Moreover, we identified the top 25 altered intercellular signaling between glial cells and neurons, highlighting their regulatory potential on gene expression in receiver cells. We reported 38 cCREs linked to sEOAD-associated genes overlapped with late-onset AD risk loci, and sEOAD cCREs enriched in neuropsychiatric disorder risk loci. This atlas helps dissect transcriptional and chromatin dynamics in sEOAD, providing a key resource for AD research.

association studies (GWAS) have identi ed over 75 common variants linked to LOAD, with the APOE ε4 allele being the most signi cant genetic risk factor [17][18][19][20][21][22] .Moreover, studies have shown that most identi ed variants associated with LOAD reside in non-coding regions and are suggested to modulate gene expression through non-direct regulatory mechanisms in a cell type-speci c manner 4,[7][8][9]23 .
Therefore, a systematic assessment of variants and their potential functions is required for understanding the genetic similarities and differences between sEOAD and LOAD.
To bridge these research gaps and delineate the molecular landscape of sEOAD, we generated singlenucleus multiome data of a small sEOAD cohort.The multiome pro ling, which combines snRNA-seq and snATAC-seq in a single experiment, can effectively reduce batch effects among experiments.Our dataset includes gene expression and chromatin accessibility pro les from over 70,000 nuclei obtained from postmortem brain samples of the PFC, EC, and HIP regions of sEOAD individuals and age-matched control donors.Here, we conducted systematic bioinformatics analyses to detect cellular compositional changes, transcriptome alterations associated with sEOAD, and cell type-speci c cCREs in uencing gene expression across different brain regions.We examined transcription factors (TFs) that may exert their functions through distal enhancer-gene connections and identi ed their targeted genes in the brain of individuals with sEOAD.Furthermore, we examined the altered intercellular communications between non-neuronal and neuronal cell types and identi ed potential target genes in receiver cells in sEOAD.
Finally, to explore shared and unique genetic components, we ne-mapped the genetic variants associated with LOAD and other traits to the cell type-speci c chromatin accessible peaks in sEOAD brains.Altogether, our study provides a comprehensive overview of gene regulatory mechanisms in sEOAD brains, shedding light on the complex interplay between diverse brain cell types in its pathophysiology.

Results
Joint snRNA-seq and snATAC-seq pro ling of sEOAD in three brain regions To investigate the dysregulated gene regulatory mechanisms in sEOAD brains, we generated joint snRNA-seq and snATAC-seq (snMultiome) pro les from each sample.We conducted 10x Chromium single-nucleus multiome assays on the nuclei isolated from postmortem PFC (n = 9), EC (n = 6), and HIP (n = 5) tissues of the sEOAD individuals and matched controls.The four individuals with sEOAD (age at death: 59-64 years) did not have any established causal variants in APP, PSEN1, and PSEN2 associated with EOAD, as documented on the Alzheimer's Forum website (https://www.alzforum.org/mutations).
The ve age-matched controls (age at death: 61-65 years) had no neuropsychiatry diagnoses (Fig. 1a & Supplementary Table 1).Sample-level quality control was performed using the CellBender tool to remove systematic background noise in the snRNA-seq assay, followed by ltration of low-quality nuclei and putative doublets in both assays (Supplementary Fig. 1, Methods).In total, we generated high-quality transcriptomic and epigenomic pro les of 71,763 nuclei.The normalization, unsupervised dimension reduction, and integration were performed based on either gene expression, open chromatin accessibility pro les, or both modalities combined.The results revealed consistent cell clustering structures across three brain regions, as visualized using Uniform Manifold Approximation and Projection (UMAP) (Fig. 1b, Supplementary Fig. 2).
We annotated the 71,763 nuclei by mapping the snRNA-seq pro le to reference snRNA-seq datasets of the human brain from previous studies 24,25 .These 71,763 nuclei were categorized into eight major cell types (with con dence score > 0.95 from references) and high-resolution sub-cell types (Fig. 1b, Supplementary Fig. 2a).Consistent with previously published snRNA-seq studies of the human brain 7,9,[24][25][26] , we observed cell type-speci c marker gene expression within each annotated cluster, including AQP4 in astrocytes, SLC17A7 in excitatory neurons, GAD2 in inhibitory neurons, CSF1R and CD74 in microglia, MOBP in oligodendrocyte, PDGFRA in oligodendrocyte precursor cells (OPCs), FLT1 in endothelial cells, and PDE5A in vascular and leptomeningeal cells (VLMC)/pericytes (Fig. 1c).Chromatin accessibility information around cell type marker genes also revealed unambiguous patterns distinguishing major cell types in the brain (Fig. 1d, Supplementary Fig. 2b-d).Together, these results indicated that the snRNA-seq and snATAC-seq data were of high quality and that cell type annotations were reliable.Because VLMC/pericytes and endothelial cells accounted for fewer than 200 nuclei, they were excluded from further analyses.
Our cell type composition analysis using the scCODA tool 27 revealed the alterations in cell type abundance among different brain regions and conditions (Fig. 1e, Extended Data Fig. 1, Supplementary Table 2).Excitatory and inhibitory neurons were signi cantly less abundant in the EC and HIP compared to the PFC (posterior probability of inclusion > 0.95), while the other six major cell types remained stable across regions (Extended Data Fig. 1a).These differences may be due to tissue sampling closer to the white matter in the HIP and EC, resulting in a higher proportion of glial cells compared to the PFC.Although a lower abundance of excitatory neurons was observed in the sEOAD PFC compared to controls, no signi cant differences in the proportion of other cell types were found between sEOAD and controls across the three brain regions (Extended Data Fig. 1b-d).These ndings aligned with the notion that AD brains are pathologically characterized by neuronal loss, at least in the PFC 28 .Although not statistically signi cant, a higher abundance of oligodendrocytes was observed in sEOAD than in controls, particularly in the EC (Extended Data Fig. 1b-d).

Cell type-speci c dysregulation of gene expression in sEOAD across brain regions
To systematically characterize transcriptomic changes, we performed differential gene expression analysis in each major cell type of different brain regions between the sEOAD and control groups.We identi ed 1,222 sEOAD-associated differentially expressed genes (DEGs) in the PFC, 3,359 in the EC, and 5,430 in the HIP (Supplementary Data 1).These sEOAD DEGs were identi ed based on consensus signals from both MAST and mixed-effect models (adjusted p-value < 0.05, absolute log2 fold change [|log2FC|] > 0.25), accounting for sex, age at death, and batch (Fig. 2a, Methods).Despite similar numbers of nuclei captured in each brain region (Supplementary Fig. 2b-d), nearly three times as many sEOAD DEGs were found in the EC and HIP compared to the PFC.This nding from the snMultiome assay supports the previous reports that the EC and HIP are primarily affected in the early stages of AD pathogenesis 29 .Furthermore, sEOAD DEGs exhibited larger effect sizes (|log2FC| >1) in the HIP and EC compared to the PFC (Fig. 2b).
Among the DEGs in the sEOAD PFC, 249 were consistent with those identi ed in previous snRNA-seq studies of LOAD PFC (Fisher's exact test, p-value < 2.2 × 10 -16 , Extended Data Fig. 2a).Notably, 184 DEGs were upregulated or downregulated in the same major cell types reported in the snMultiome study by Anderson et al. 9 .Three DEGs were consistent across all studies: upregulation of MRAS in astrocytes and upregulation of CHORDC1 and TP53TG5 in oligodendrocytes.
We identi ed 362 cell type-speci c sEOAD DEGs consistently upregulated and downregulated across all three brain regions (Fig. 2c, Extended Data Fig. 2b), including 88 genes previously reported in snRNA-seq studies of LOAD PFC 7,9,24 .For example, CACNA1A was consistently upregulated in microglia across the PFC, EC and HIP (Fig. 2d, e, Supplementary Fig. 3).CACNA1A encodes a subunit of a calcium-dependent voltage channel (CaV2.1).The upregulation of this gene, speci cally in microglia, was reported to be associated with Aβ load in the human brain 30 .We also identi ed other consistently dysregulated genes, such as downregulated CIRBP and upregulated ITGB8 in the astrocytes (Extended Data Fig. 2c).Those genes have been reported to be relevant to AD mechanisms 31,32 .
Among the consistent DEGs across regions, we reported 274 genes that were not previously reported in snRNA-seq analysis of LOAD in the PFC.For example, CXXC4, which we found upregulated in excitatory neurons, is part of the Wnt signaling pathway and a key regulator of the amyloid precursor protein (Extended Data Fig. 2c).MAML2, which we found downregulated in oligodendrocytes, is part of the Notch signaling pathway, linked to neurovascular dysfunction in AD 33 .Additionally, 86 sEOAD DEGs were consistently downregulated in inhibitory neurons across three brain regions.Gene set enrichment analysis (GSEA) revealed their enriched Gene Ontology (GO) Biological Process (BP) terms, such as "modulation of chemical synaptic transmission", "regulation of trans-synaptic signaling", and "synapse organization" [Benjamini-Hochberg (BH) adjusted p-value < 0.05, Extended Data Fig. 2d].
Most cell type-speci c sEOAD DEGs showed varying dysregulation patterns across the considered regions.For example, GRM3, which encodes the glutamate metabotropic receptor 3 and is associated with synaptic function, was downregulated in astrocytes in the PFC and EC but not in HIP (Fig. 2f).This gene was found altered in AD and other neurological disorders 34,35 .
In our GSEA of sEOAD DEGs for each cell type, we observed diverse enriched GO BP terms across the three regions (BH adjusted p-value < 0.05) (Fig. 2g).Upregulated DEGs in astrocytes and inhibitory neurons were consistently enriched in synaptic functions, such as "modulation of chemical synaptic transmission" and "regulation of trans-synaptic signaling," across all brain regions.However, only upregulated DEGs in excitatory neurons in the HIP were enriched in synaptic functions, indicating regionspeci c dysregulation.Upregulated genes in the microglia were enriched in cell morphology-related functions in the EC and HIP but not in the PFC.Downregulated DEGs in inhibitory neurons across all brain regions and in excitatory neurons in the HIP were enriched in "oxidative phosphorylation," "aerobic respiration," and "ATP metabolic process," suggesting potential mitochondrial dysfunction and reduced energy metabolism in sEOAD brains 36,37 .

Cell type-speci c candidate cis-regulatory elements (cCREs) for sEOAD DEGs
To assess potential cell type-speci c open chromatin accessible signals associated with sEOAD, we performed peak calling on snATAC-seq assay using MACS2.We identi ed a total of 316,172 peaks in all major cell types, 75.3% of which were overlap with the results from a previously published cCRE dataset from the normal human brain by Li et al. 4 , and 61.5% of the peaks overlapped with the signals in the ENCODE cCREs database (Fig. 3a) 38 .Among these peaks, 22.45% were within 3 kbp of the nearest transcriptional start sites (TSSs), 46.21% were in intronic regions, and 24.13% were in distal intergenic regions (Fig. 3b).
To decode the dysregulated transcriptional regulatory circuitry in sEOAD brains, we conducted peak-togene linkage analyses in each brain region.We leveraged shared barcodes in snMultiome pro ling to directly link epigenomic pro les with the respective transcriptomic pro les.We prioritized open chromatin peaks within 500 kbp of the TSS of each sEOAD DEG, focusing on signi cant positively correlated peaks (Spearman's correlation > 0.05, BH adjusted p-value < 0.05) with the target gene, considering peak size, GC content, and fragment count.In total, we identi ed 3,794 cCRE-DEG associations in the PFC, involving 724 unique DEGs.In the EC, we found 21,349 of these associations, impacting 2,437 unique DEGs; while in the HIP, we reported 27,666 associations affecting 3,171 unique DEGs.Each cCRE was linked to a median of four DEGs in the PFC, six in the EC, and ve in the HIP (Fig. 3c).The median distance between the cCRE and the TSS of the linked DEG was 140,980 bp in all three brain regions.The greatest proportion of the cCREs linked to sEOAD DEGs were present within distal enhancer-like regions based on the annotations of ENCODE cCRE database (52.00% in PFC, 51.99% in EC, and 55.11% in HIP) (Fig. 3d, Supplementary Data 2).This pattern was consistent in cCRE linked to sEOAD DEGs across cell types and brain regions (Extended Data Fig. 3).The results suggested the critical role of distal enhancers in gene expression alterations in sEOAD brains.
Transcription factors are the primary regulators of gene expression, and their dysregulation has been associated with AD pathogenesis 39 .Using snATAC-seq data, we conducted motif enrichment analysis to nominate TFs that may contribute to altered gene expression in sEOAD in each cell type of each brain region.To prioritize bona de signals, we focused on TFs expressed in over 25% of the corresponding cell types.In total, 217 TF motifs were enriched in the cCREs linked to sEOAD DEGs (BH adjusted p-value < 0.05, Fig. 3e, Supplementary Data 3).Among them, 148 motifs were signi cantly enriched in the cCREs linked to sEOAD DEGs across multiple cell types.For instance, motifs of ZNF148 had signi cant enrichment in cCREs associated with sEOAD DEGs in all major cell types (BH adjusted p-value < 0.05, Fig. 3e).Other TFs, such as EGR1/3 in excitatory neurons and SPI1 and IKZF1 in microglia, appeared to exert their regulatory functions in speci c cell types.The motif activity of these TFs within speci c cell types and brain regions was further supported by DNA footprinting analyses (Fig. 3f, g).Importantly, most identi ed TF motifs (198 out of 217) were enriched in the cCREs positively correlated with the upregulated DEGs, and they also enriched in distinct cCREs correlated with downregulated DEGs in the same cell type.This suggested TFs may function as both activators and repressors in different chromatin regions, thereby targeting different genes.Indeed, it has been reported that the effect of TFs on transcription is context-dependent; TFs can recruit different coactivators or corepressors, forming multi-subunit protein complexes to regulate transcription 40,41 .

Cell type-speci c regulatory networks in the sEOAD brain
To further explore the alterations in sEOAD brains, we identi ed cell type-speci c TFs potentially modulating gene expression through enhancer binding in each brain region (Fig. 4a).Using SCENIC + algorithm, we inferred enhancer-driven regulons for each brain region, which comprising candidate TFs, open chromatin regions with enriched TF-binding motifs, and target genes across cellular states 42 .We prioritized the regulons exhibiting the highest 5% positive correlations and the lowest 5% negative correlations between TF expression and open chromatin region-based regulon enrichment scores.We identi ed 16, 31, and 17 key regulators in the PFC, EC, and HIP, respectively (Fig. 4b, Extended Data Fig. 4a, b, Supplementary Data 4).These regulators may act as either activators or repressors.On average, each regulon comprised 1,014 open chromatin regions and 262 target genes.We observed high cooperativity in prioritized regulons, as evidenced by the shared target regions and genes among them (Extended Data Fig. 4c-e).For instance, the chromatin regions and genes targeted by IKZF1 and FLI1 were highly overlapped in all three brain regions (Jaccard coe cient > 0.5, Extended Data Fig. 4c-e).These results suggested that a complex regulatory network, involving cooperativity from multiple TFs, is necessary to in uence the transcriptome pro les in the human brain 43,44 .Seven regulators were found to be conserved in speci c cell types across all three brain regions, including RFX4 in astrocytes and FLI1, IKZF1, RUNX1, FOXN3, ETS2, and POU2F2 in microglia (Extended Data Fig. 4f).
We further explored the gene regulatory networks in the EC region.SCENIC + identi ed 31 regulons, indicating high regulatory activity (Fig. 4b).Among these regulons, we prioritized TFs with the motifs signi cantly enriched in cCREs linked to sEOAD DEGs in the corresponding cell types (Fig. 3e, Supplementary Data 3).This process pinpointed ve activators in astrocytes (RFX4, PAX6, RFX2, EMX2, and TCF7L2) and two in microglia for constructing gene regulatory networks.The ve activators in astrocytes targeted a total of 102 upregulated and 113 downregulated DEGs in the EC (Fig. 4c, Extended Data Fig. 4e).Particularly, RFX4 showed high speci city and regulatory activity in astrocytes across sEOAD and control samples, as validated by DNA footprinting analysis (Fig. 4d).The inferred chromVAR motif activity of RFX4 was signi cantly higher in astrocytes of sEOAD individuals (BH adjusted p-values = 9.17 × 10 − 32 , Fig. 4d), indicating altered regulatory activity 45 .The RFX4 regulon contained 661 target genes, including 82 genes upregulated and 96 genes downregulated in astrocytes of sEOAD individuals.RFX4 (Regulatory Factor X 4) may play a role in psychiatric conditions in regulating genes crucial for brain development and function, including those involved in circadian rhythms and neural differentiation 46 .However, the disrupted regulatory effect of RFX4 on the transcriptome alterations in astrocytes of sEOAD has not been fully explored.In our sEOAD data, we observed the upregulation of GFAP, which encodes for glial brillary acidic protein, which is associated with astrocyte activation during stress and injury 47 .The GSEA of the DEGs targeted by RFX4 in astrocytes revealed diverse enriched GO BP terms, including "L-aspartate transmembrane transport", "regulation of membrane potential", and "Wnt signaling pathway", among others (BH adjusted p-values < 0.05, Fig. 4e, f).This signi ed its potential role in regulating multiple transport-related pathways in the brain of sEOAD.Furthermore, we prioritized activators in the microglia of EC, including FLI1 and IKZF1.Both TFs exhibited high speci city within microglial cells, as their motifs were signi cantly enriched in cCREs linked to sEOAD DEGs in microglia.Speci cally, these TFs targeted 71 upregulated and 20 downregulated DEGs (Fig. 4g).DNA footprinting analysis of IKZF1 showed high motif activity in the microglia of both sEOAD and control samples.Moreover, the inferred chromVAR motif activity was notably increased in sEOAD (adjusted p-value = 0.03, Fig. 4h).The IKZF1 regulon involved 47 upregulated and 20 downregulated DEGs in the microglia of sEOAD.The dysregulated genes targeted by IKZF1 were predominantly involved in immune-related functions such as "negative regulation of immune response", "leukocyte-mediated immunity", and "immune response-activating signaling pathways" (BH adjusted pvalue < 0.05, Fig. 4i).This underscores the critical role of IKZF1 in mediating the immune response and pathophysiology of sEOAD.
Collectively, these results suggested that the cell type-speci c transcriptomic pro les were regulated through cooperative interactions of multiple TFs and enhancers on their target genes in sEOAD brains.We highlighted a disrupted regulatory role for RFX4 in the astrocytes of sEOAD EC, which modulates genes involved in intercellular signaling pathways, and for IKZF1 in microglia, related to innate immunity and neuroin ammation.

Dysregulated intercellular signaling between non-neuronal and neuronal cell types
The altered regulatory relationships in non-neuronal cell types and their impact on neuronal circuit homeostasis and pro-in ammatory responses might play an important role in LOAD brains 48 .To estimate such association in sEOAD brains, we conducted comparative cell-cell communication analyses using the MultiNicheNet tool 49 .This analysis facilitated the identi cation of dysregulated intercellular communication mediated by ligand-receptor pairs and their potential regulatory in uence on the receiver cells.We prioritized the top 25 differential signals sent and received by speci c cell types as condition-speci c intercellular signaling.
Here, we highlighted the intercellular communication analysis of astrocytes and microglia in the EC.Speci cally, the results showed that all top 25 differential signals sent by astrocytes were downregulated in sEOAD (Fig. 5a).The communication signals between EFNA5 in astrocytes to EPHA3 and EPHA6 in excitatory neurons were downregulated (Fig. 5a).The Ephrin-Eph signaling plays an important role in regulating cell migration in neurons.It is also associated with cytoskeletal rearrangements, as it regulates small GTPase activation and inactivation 50 .Furthermore, we observed multiple downregulated ligand-receptor pairs involved in neuroin ammation regulation, such as PTN-PTPRS/PTPRB signaling from astrocytes to excitatory neurons and A2M-LRP1 signaling from microglia to astrocytes (Fig. 5a, b).Remarkably, although more upregulated intercellular signalings were observed in the sEOAD PFC, downregulated PSAP-LRP1 signaling from microglia to astrocytes was noted in all three brain regions (Extended Data Fig. 5a-d).The reduced levels of PSAP protein have been associated with increased neuro brillary tangle development and neuroin ammation 51 .Furthermore, we identi ed upregulated TGFB1-ITGB8 signaling from microglia to astrocytes, suggesting a potential increasing microgliogenesis in the sEOAD EC and PFC 52 .
We further explored the regulatory in uence of ligands in sender cell types on the sEOAD DEGs in receiver cell types.Among the top differential communication pathways from astrocytes to excitatory neurons, ligands encoded by the EFNA5, OMG, and PTN in astrocytes exhibited increased ligand activity (Fig. 5c).Ligand EFNA5 exhibited the highest regulatory potential on downregulated expression of ARC.The Activity-Regulated Cytoskeleton-associated protein encoded by ARC plays a signi cant role in synaptic plasticity and neuronal activity 53 .Additionally, the TF coding genes in excitatory neurons, such as EGR1, EGR3 and FOS, were targeted by ligands in astrocytes (Fig. 5c).Similarly, the ligands EFNA5 and PTN in astrocytes displayed elevated regulatory activities on downstream genes in excitatory neurons in the PFC (Extended Data Fig. 5e).Conversely, ligands encoded by PSAP and A2M in microglia exhibited high ligand activity, potentially modulating the expression of targeted genes in astrocytes (Fig. 5d, Extended Data Fig. 5e, f).These results underscored the potential subtle and dynamic intercellular signaling in uences on neuronal and glial functions in sEOAD brains.

Genetic variation at candidate CREs linked to sEOAD DEGs
To better understand the cellular impact of genetic variants identi ed by AD GWAS, we prioritized risk loci that overlapped with cCREs linked to sEOAD DEGs.We curated 148 index single nucleotide polymorphisms (SNPs), reported as genome-wide signi cant from four LOAD GWA studies (Fig. 6a) 17,18,20,21 .We expanded these SNPs using linkage disequilibrium (LD) to include nearby variants with high coinheritance probability (LD R 2 > 0.8 based on phase 3 of the 1000 Genomes Project data).We identi ed 3,180 unique SNPs associated with LOAD, of which 98.2% were in non-coding regions.Interestingly, 133 LOAD SNPs mapped in cCRE sequences were linked to 29 sEOAD DEGs across each brain region (Extended Data Fig. 6a).Particularly, cCREs in microglia in the EC and HIP were signi cantly enriched with LOAD GWAS variants, unlike those in the PFC (adjusted Fisher's exact test, p-value < 0.05, Fig. 6b).In addition, cCREs linked to sEOAD DEGs of excitatory neurons in HIP were enriched with LOAD GWAS loci.Notably, cCREs linked to the downregulated HLA-DRB1 and HLA-DRA genes in microglia intersected with the largest number of SNPs associated with LOAD (Extended Data Fig. 6a).The index SNP rs9271058 (reported in Kunkle et al. 18 ) was mapped to a cCRE (chr6:32607513-32607728) linked to HLA-DRB1 in microglia in the EC (Fig. 6c).There were 43 SNPs mapped to a cCRE (chr6:32603881-32604799) that was linked to HLA-DRB1.Importantly, this distal-like enhancer cCRE contains TF motifs speci c to microglia regulators, such as IKZF1 and FLI1.We found that the index SNP (rs74685827), and its associated SNPs reported by Bellenguez et al. 21, intersected with the microglia-and astrocytes-speci c cCREs linked to dysregulated SORL1 gene expression in the HIP (Extended Data Fig. 6b).We extended this approach to include 172 unique SNPs (LD R 2 > 0.8) from ve index SNPs associated with EOAD (p-value < 5×10 − 8 ) 54 .We pinpointed 20 SNPs that intersected with cCREs correlated with HLA-DRB1 and MS4A6A in microglia within the EC, and 16 SNPs intersected with cCREs linked to HLA-DRB1 and HLA-DRA in the HIP.We further examined whether the cCREs of sEOAD DEGs harbored single-cell expression quantitative trait locus (sc-eQTL) identi ed in the human brain 55 .We identi ed 6, 27, and 33 brain sc-eQTL loci intersecting with cCREs of sEOAD DEGs in the PFC, EC, and HIP, respectively (Extended Data Fig. 6c, d).
Similarly, we performed sLDSC analysis on cell type-speci c cCREs (Methods).As expected, microgliaspeci c cCREs demonstrated signi cant enrichment for the heritability of LOAD, as well as various immune-mediated disorders, such as Crohn's disease, multiple sclerosis, rheumatoid arthritis, ulcerative colitis, and vitiligo (adjusted p-value < 0.05, Extended Data Fig. 6e, Supplementary Data 5b).A recent study revealed shared and unique genetic components among these disorders 75 .Furthermore, cell typespeci c cCREs in excitatory neurons were signi cantly enriched in risk variants associated with neuropsychiatric disorders.
In summary, the cell type-speci c cCREs linked to sEOAD DEGs suggested that LOAD risk variants in open chromatin peaks may account for a limited number of microglia transcriptome alterations.Further research is needed to determine the distinct genetic risk factors related to sEOAD.

Discussion
We reported the rst single-nucleus multiome pro ling of sEOAD, creating a comprehensive atlas of gene expression and regulation in the PFC, EC, and HIP regions.By pro ling 71,163 nuclei from four sEOAD cases and ve matched controls, we simultaneously assessed chromatin accessibility and gene expression.Our analyses identi ed cell type-speci c transcriptome alterations in sEOAD and linked cCREs across brain regions.We characterized conserved regulons in non-neuronal cells and found that multiple TFs may cooperatively modulate transcriptomic changes, leading to dysregulated synaptic plasticity and neuroin ammation in astrocytes and microglia.Additionally, we prioritized altered intercellular signaling pathways and their regulatory impacts on sEOAD DEGs.Our genetic ne-mapping analysis revealed cCREs linked to sEOAD-associated genes intersected with late-onset AD risk loci and may share genetic risks with neuropsychiatric disorders.
The inclusion of multiple brain regions in our study has enabled us to examine conserved regulators and their target genes across the brain regions that are disrupted in sEOAD.Our results highlighted regulators with high speci city in glial cells, such as RFX4 in astrocytes and IKZF1 in microglia.Remarkably, the RFX4 motif was enriched in cCREs linked to sEOAD DEGs across all regions, with targeted DEGs consistently enriched in transmembrane transport and trans-synaptic signaling functions.This suggests a conserved role for RFX4 in astrocytes and its universal regulatory alterations in sEOAD pathogenesis.Similarly, in microglia, upregulated genes related to neuroin ammation were targeted by TFs like IKZF1 and FLI1 in sEOAD brains.Additionally, we found cCREs contain TF motifs speci c to microglia regulators, such as IKZF1 and FLI1, suggesting that SNPs in these cCREs may impact gene expression regulation to confer sEOAD risk.Our ndings indicate that gene expression changes in sEOAD are regulated by complex synergistic and/or antagonistic effects of multiple TFs through distal enhancers.
We further prioritized cCREs linked to sEOAD DEGs in each brain region through integrated epigenomic and transcriptomic analysis.Interestingly, we found that cCREs linked to sEOAD DEGs in microglia intersected with genetic risk loci identi ed in LOAD GWAS studies 17,18,20,21 .However, we did not nd signi cant heritability enrichment in the sLDSC analysis (Fig. 6d).Genetic variants associated with psychiatric disorders and traits, such as schizophrenia, bipolar disorder, and neuroticism, were signi cantly enriched in both the cCREs linked to sEOAD DEGs and cell type-speci c cCREs in excitatory neurons.This suggested potential shared genetic risks between sEOAD and psychiatric conditions, particularly in excitatory and inhibitory neurons 4,7 .Our integrated epigenomic and transcriptomic analysis successfully explored the shared and unique genetic components and their potential functions associated with sEOAD.
Our study offers insights into cell type-and brain region-speci c transcriptome dysregulation and potential regulatory disturbances in sEOAD.However, our study has limitations.First, this study is constrained by a small sample size, resulting from the limited availability of postmortem brain samples in sEOAD.We conducted a simulation-based power analysis and found lower statistical power estimated for identifying differentially expressed genes in less abundant cell types, such as inhibitory neurons (Supplementary Table 4).This limitation, due to the availability of brain tissue and also the high cost in multiome, also restricts our ability to perform further analyses at the sub-cell type level.We anticipate that larger-scale single nuclei studies will enhance the accuracy and detail in identifying gene expression changes in sEOAD.Second, future laboratory experiments should help validate the candidate TFs and cCREs in speci c cell types.Third, our study did not comprehensively assess the regulatory mechanisms involving the complex interactions of multiple regulators on target genes.Finally, we lacked the access to sEOAD GWAS datasets for combined association with our scMultiome features, since they are not publicly available yet.We will integrate such analysis when sEOAD GWAS datasets are released to the public.
In conclusion, our study provides novel insights into transcriptome alterations and identi es key gene regulators, such as RFX4 in astrocytes and IKZF1 in microglia, that are crucial in the pathophysiology of sEOAD across three brain regions commonly affected in sEOAD.The intercellular communication analysis revealed the altered roles of non-neuronal cell types in synaptic plasticity, neuroin ammation, and microgliogenesis in sEOAD.We delineated distinct genetic architectures within cell type-speci c cCREs.These discoveries set the stage for future investigations into the unique mechanisms of sEOAD pathophysiology and for the development of targeted therapeutic strategies.

Human postmortem brain sample collection
Postmortem brain samples of nine individuals were obtained from the NIH NeuroBioBank from the University of Miami and McGovern Medical School at UTHealth (Supplementary Table 1) in compliance with UTHealth Houston's Institutional Review Board.sEOAD cases from the McGovern Medical School were de ned as subjects diagnosed with sEOAD before the age of 65 years and con rmed neuropathological diagnosis postmortem.Brain samples of sEOAD from the NIH NeuroBioBank were selected based on the neuropathological diagnosis of sEOAD with onset before age 65 years.The genotype of the samples from the NIH NeuroBioBank were examined based on paired whole-genome sequencing data publicly available on The National Institute of Mental Health Data Archive (https://nda.nih.gov/edit_collection.html?id=3917).Subjects with known autosomal dominant mutations in APP, PSEN1, and PSEN2 were excluded from the analyses.The sEOAD cases from McGovern Medical School were the subjects diagnosed with sEOAD before the age of 65 years and con rmed neuropathological diagnosis postmortem.All tissues were de-identi ed under the privacy rules of the Health Insurance Portability and Accountability Act of 1996 (HIPAA).The tissue donation consent was obtained from all participants by UTHealth Institutional Review Board.

Tissue preparation and nuclei isolation
The tissue preparation, dissociation, and nuclei extraction were performed at the UTHealth Cancer Genomic Core.Samples were stored at -80°C.The nuclei isolation protocol was modi ed based on a published protocol 76 and 10x Genomics protocol (CG000375).Brie y, 30 mg of frozen tissue was nely minced on dry ice and then homogenized in buffer on ice within 20 minutes.The homogenized tissue was ltered through a 35 µm cell strainer.Small cell debris was removed by centrifuging at 300g for 3 minutes at 4°C.The large cell debris was removed by gradient centrifugation using OptiPrep density gradient medium (Sigma, cat# D1556-250ML) at 7200g for 15 minutes at 4°C.The nuclei pellet was resuspended in a 0.1x lysis buffer.After washing the isolated nuclei in the wash buffer, the nuclei were resuspended in 1x nuclei buffer before single-nuclei Multiome library construction.

Single-nuclei multiome sequencing
Single-nuclei libraries were constructed following the 10x Chromium Single Cell Multiome ATAC + Gene Expression protocol (CG000338).Brie y, nuclei suspensions were incubated with a transposase, which fragmented the DNA in open regions of the chromatin and added the adapter sequences to the ends of the DNA fragments.The transposed nuclei were loaded onto Chromium Next GEM Chip J (PN-1000234, 10x Genomics, Pleasanton, CA) with portioning oil and barcoded single-cell gel beads, followed by PCR ampli cation.The ATAC and the gene expression libraries were then prepared separately.Library qualities were assessed using the Agilent High Sensitive DNA Kit (#5067 − 4626) on an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, USA).Quali ed libraries underwent paired-end sequencing on an Illumina NovaSeq System (Illumina, Inc., USA).

Joint single-nuclei multiome data processing and cell type annotation
We utilized the Cell Ranger ARC (v2.0.2) pipeline for aligning the snRNA-seq and snATAC-seq reads from each sample to the GRCh38 genome and quantifying them with 'cellranger-arc count'.The raw snRNA count matrix and snATAC fragment counts were further processed using a rigorous pipeline.Firstly, we removed the systematic background noise in the raw snRNA count matrix using CellBender, which applies a deep-learning model to accurately distinguish cell-containing droplets from cell-free droplets 77 .Secondly, we applied initial ltration based on RNA-assay matrices (200 < nFeature_RNA < 7,500) and ATAC-assay matrices (Nucleosome_signal < 2.5, TSS.enrichment > 2) using Seurat (v5.0.0) 78 and Signac (v1.12.0) 79 packages, respectively.Potential doublets were removed using the DoubletFinder R package using default parameters 80 .Then, the RNA counts were normalized using SCTransform with mitochondrial percent per cell regressed out.Principal component (PC) analysis (PCA) and UMAP visualization based on the rst 30 PCs were performed on normalized RNA counts.Additionally, we performed pseudo-bulk analysis on RNA assay to identify and remove outlier samples through PCA.
Predicted major cell type annotations were determined for each cell based on SCTransform normalized reference mapping embedded in Seurat.Gene expression (RNA assay) was mapped to SCTransformed normalized snRNA-seq datasets of the PFC and middle temporal gyrus 24,25 and annotated with the cell types of the reference map.Cells with consistent cell type predictions and predicted cell type scores greater than 0.95 in both datasets were kept for downstream analysis.Additionally, we leveraged scANVI from scvi-tools (v1.0.3) for high-resolution cellular subclasses annotation based on the reference snRNA-seq dataset 25,81 .
After predicting the major cell types based on RNA assay, we performed open chromatin peak calling on ATAC assay at pseudo-bulk level with MACS2 82 using the CallPeaks function in Signac.Peaks that overlapped with genomic blacklist regions for the GRCh38 genome were removed, and peak counts were normalized using Latent Semantic Indexing (LSI), including term-frequency inverse-document-frequency (TFIDF).Dimension reductions were performed with singular value decomposition (SVD) of the normalized open chromatin accessibility matrix.The open chromatin accessibility UMAP visualization was created using the second through the 50th LSI components.
To mitigate batch effects, we applied the Harmony (v1.0.0) 83 algorithm on both snRNA-seq and snATACseq data.Then, we created a weighted nearest neighbor (WNN) graph with the FindMultiModalNeighbors function of Seurat to represent a weighted combination of both modalities.

Cell type proportional analysis by scCODA
We applied the single-cell compositional data analysis using scCODA (v0.1.9) 27, a Bayesian model, to detect the cell type proportional shift across brain regions and conditions.The statistical test and the posterior inclusion probabilities (Pinc) for credible effects were calculated by including either region or sEOAD diagnosis labels as covariates at a false discovery rate (FDR) < 0.05.

Differential expression analysis
The DEGs between sEOAD and control samples (sEOAD DEGs) were determined for each major cell type across brain regions.VLMC/pericytes and endothelial cells were excluded from the DEG analysis due to low cell counts.First, we conducted DEG analysis using the MAST framework for the genes presented in at least 25% nuclei in either of the two condition groups.Gene expression data underwent lognormalization with a scale factor of 1 × 10 5 .Age, sex, and batch were included in the MAST model as covariates.Second, we implemented a mixed-effect model, utilizing the Libra R package 84 , to account for the random effect of subject origin, which accommodates a zero-in ated negative binomial distribution.The proposed formula of the mixed-effect model is: Diagnosis + Sex + Age at death + Batch + offset(log(total counts)) + (1|subject) Genes were considered sEOAD DEGs if they met a Bonferroni adjusted p-value threshold of < 0.05 in both the MAST and the mixed-effect model, along with an absolute log2FC greater than 0.25 between sEOAD and control in speci c cell types.

Gene Ontology (GO) enrichment analysis
ClusterPro ler (v4.10.0) was used for gene set enrichment analysis to nd potential GO BP terms for upor down-regulated sEOAD DEGs in each cell type within each brain region 85 .The background gene set for this analysis included all genes listed in the GO database.Additionally, we used the simplify function to remove redundant enriched GO BP terms based on calculated similarity (similarity score > 0.8).The multiple testing correction was performed, and FDR < 0.05 was applied.

Open chromatin peaks annotations
The open chromatin peaks identi ed in each major cell type were annotated using ChIPseeker 86 and TxDb.Hsapiens.UCSC.hg38.knownGeneusing default 3-kb upstream and 3-kb downstream of the TSS.

GeneExpression ∼
Additionally, we examined the intersection between open chromatin peaks identi ed in our study and three published datasets: cCRE identi ed in single nuclei from the adult human brain 4 and ENCODE Registry of cCRE in the human genome 38 .Brain cell type-speci c cis-eQTL data was obtained from Bryois, J. et al. 55

Peak-gene linkage analysis
We conducted peak-gene linkage analysis on sEOAD DEGs resulting from the previous step in each cell type using the LinkPeaks function in Signac 79 , based on the approach initially described by SHARE-seq 87 .The Spearman correlation coe cient between gene expression and peak accessibility was calculated.
Open chromatin peaks within 5 × 10 5 bp from the gene TSS were included in the model.The GC content, overall accessibility, and peak size were considered in the model as covariates for correcting bias.Benjamini-Hochberg (BH) adjusted p-values were calculated in each region.Only high-con dence positive peak-gene links with BH-adjusted p-value < 0.05 and correlation coe cients > 0.05 were retained for downstream analyses.

Single-nuclei TF binding motif analysis
We retrieved a comprehensive of TF binding motif position weight matrices covering 755 TFs from the JASPAR 2020 database (collection = "CORE", species = "Homo sapiens", access date: March 9th, 2024).The DNA sequence motif information was added to the snATAC-seq assay using the AddMotifs function, and the motif activity matrix was calculated using chromVAR implemented in the Signac package 79 .
To prioritize the TFs associated with sEOAD, we identi ed overrepresented motifs for each set of cCREs linked to sEOAD DEGs in different cell types and brain regions.This involves employing the FindMotifs function to conduct a hypergeometric test, evaluating the likelihood of the observed motif frequency occurring by chance relative to a GC content-matched background peak set.The background set was uniquely determined for each set of peaks linked to DEGs in each cell type.Enriched TFs were kept if the TF gene was expressed in over 25% of the corresponding cell type.
Additionally, we identi ed TFs with differential motif activity by performing differential testing on the chromVAR z-score between sEOAD and controls in each major cell type using the FindMarker function at Bonferroni adjusted p-value < 0.05.Footprinting analysis of selected TFs was conducted, accounting for the Tn5 sequencing insertion bias.

Gene regulatory network inference
We applied the SCENIC + work ow to construct TF regulatory networks in speci c cell types and conditions 42 .Firstly, the topic modeling of the snATAC-seq data was carried out using pycisTopic to identify the sets of co-accessible regions and the candidate enhancer regions.Next, the motif enrichment analysis on the candidate enhancer regions was carried out using cis-target and DEM algorithms on the cell type-and condition-speci c differentially accessible regions.Lastly, the TF regulatory networks were built by calculating the accessible region-gene and TF-gene relationships using a gradient-boosting machine regression method.The retrieved regulatory networks were ltered out if the region-to-gene relationship was negative.The TF-cistrome relationships (Rho) were calculated on the ltered TF regulatory networks, and only the regulons with the top 5% highest positive correlation values and 5% lowest negative correlation values were retained for further downstream analysis.

Differential cell-cell communication analysis
We employed MultiNicheNet (v1.0.1, last updated on March 20th, 49 , a universal tool to infer the putative context-speci c cell-cell communication signals mediated by ligand-receptor pairs across cell type and brain regions in sEOAD and control samples.MultiNicheNet ranks the importance of ligandreceptor pairs based on the 'pseudo-bulk' transcriptome pro le of each sample.Potential ligands were extracted if expressed in at least 25% of the sender cells within their respective clusters, and the p-value cutoff was set to 0.05.To prioritize highly variable intercellular signals in sEOAD, we highlighted the top 25 differential signals sent and received by selected cell types per experimental group.Furthermore, we applied the make_ligand_activity_target_plot function and calculated the regulatory potential scores between prioritized ligands and potential targeted genes within the receiver cell using default parameters.

Genetic variants selection and LD expansion
Index SNPs were collected from four large-scale LOAD GWA studies 17,18,20,21 .LD expansion was performed using LDlinkR 88 to identify SNPs with LD R 2 > 0.8 based on the Phase 3 of The 1000 Genomes European ancestry 89 .We performed Fisher's exact test to assess the enrichment of genetic variants associated with LOAD in cCREs linked to sEOAD DEGs in each cell type across brain regions with a Holm-Bonferroni adjusted p-value < 0.05 as the cutoff.The background cCREs set was de ned as open chromatin peaks identi ed in each cell type but not linked to sEOAD DEGs.
We performed sLDSC in two scenarios.First, the sLDSC analysis was conducted on cCREs linked to sEOAD DEGs in each cell type in each region (Spearman correlation > 0.05; BH adjusted pvalues < 0.05).Second, we applied sLDSC on cell type-speci c cCREs.The cell type-speci c cCREs were identi ed using the FindAllMarkers function from the Seurat package.This function employed a logistic regression model and adjusted for the library size.We prioritized the peaks exhibiting an Bonferroni adjusted p-value < 0.01 and absolute log2FC ≥ 1 in each cell type as the cell type-speci c cCREs.

Power analysis of differentially expressed genes analysis
To estimate the power of our dataset in detecting alterations in gene expression, we performed a simulation-based power analysis using the R package Hierarchicell 90 .This algorithm takes multiple factors into consideration, including gene dropout rates, intra-individual dispersion, inter-individual variation, variable or xed number of cells per individual, and the correlation between cells within an individual.We employed compute_data_summaries and power_hierarchicell functions to simulate and estimate the power with the following parameters: the number of samples in each condition (n = 5, PFC; n = 3, EC; n = 3, HIP), the minimal fold change between conditions (rho = 1.5 or 2, two parameters), the pvalue cutoff (p = 0.05), mean of the number of cells in each sample of each condition, and number of genes (n = 3000).

Statistics and
All statistical and analyses used in this study were described in the Methods, gure legends, or main text as appropriate.

Figure 1 Single
Figure 1